In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
In [9]:
xlow = 0
xhigh = np.pi/2
x = np.linspace(xlow, xhigh, num = 100)
y = np.cos(x)
In [10]:
plt.figure()
plt.plot(x, y)
plt.show()
In [18]:
dim = 100
x_mcmc = np.random.random(dim) * (xhigh - xlow) + xlow
y_mcmc = np.cos(x_mcmc)
In [19]:
plt.figure()
plt.plot(x_mcmc, y_mcmc, 'o')
plt.show()
In [21]:
y_mcmc.mean() * (xhigh - xlow)
Out[21]:
In [ ]:
(xhigh - xlow) * np.sqrt(y_mcmc.
In [1]:
from numpy.random import random
the following code is from http://code.activestate.com/recipes/414200-metropolis-hastings-sampler/
In [23]:
def sdnorm(z):
"""
Standard normal pdf (Probability Density Function)
"""
return np.exp(-z*z/2.)/np.sqrt(2*np.pi)
n = 10000
alpha = 1
x = 0.
vec = []
vec.append(x)
innov = random(n) * 2 * alpha - alpha #random inovation, uniform proposal distribution
for i in xrange(1,n):
can = x + innov[i] #candidate
aprob = np.min([1.,sdnorm(can)/sdnorm(x)]) #acceptance probability
u = random(1)
if u[0] < aprob:
x = can
vec.append(x)
In [24]:
len(vec)
Out[24]:
In [25]:
#plotting the results:
#theoretical curve
x = np.arange(-3,3,.1)
y = sdnorm(x)
plt.subplot(211)
plt.title('Metropolis-Hastings')
plt.plot(vec)
plt.subplot(212)
plt.hist(vec, bins=30,normed=1)
plt.plot(x,y,'ro')
plt.ylabel('Frequency')
plt.xlabel('x')
plt.legend(('PDF','Samples'))
plt.show()
my version of the code
In [41]:
def metropolis_hastings(f, init = 0, nsamples = 1000):
alpha = 1
sample_count = 0
x = init
vec = []
vec.append(x)
while sample_count < nsamples-1:
can = x + random(1)[0] * 2 * alpha - alpha #candidate
aprob = np.min([1.,sdnorm(can)/sdnorm(x)]) #acceptance probability
u = random(1)
if u[0] < aprob:
x = can
vec.append(x)
sample_count+=1
return np.array(vec)
In [42]:
vec = metropolis_hastings(sdnorm, nsamples=10000)
In [43]:
len(vec)
Out[43]:
In [44]:
#plotting the results:
#theoretical curve
x = np.arange(-3,3,.1)
y = sdnorm(x)
plt.subplot(211)
plt.title('Metropolis-Hastings')
plt.plot(vec)
plt.subplot(212)
plt.hist(vec, bins=30,normed=1)
plt.plot(x,y,'ro')
plt.ylabel('Frequency')
plt.xlabel('x')
plt.legend(('PDF','Samples'))
plt.show()
In [45]:
def gauss(z):
"""
Standard normal gaussian
"""
return np.exp(-z*z/2.)
In [48]:
import scipy.integrate as integrate
In [50]:
integrate.quad(gauss, -np.Inf, np.Inf)
Out[50]:
In [51]:
np.sqrt(2*np.pi)
Out[51]:
In [167]:
def triangle(z):
if (type(z) is not type([1])) and (type(z) is not type(np.array(1))):
x = np.array([z])
else:
x = np.array(z)
return np.piecewise(x, [x < 0, (x < 10) & (x > 0), (x >= 10) & (x <= 20), x > 20], [0, lambda x: x, lambda x: - x + 20, 0])
In [168]:
triangle([1])
Out[168]:
In [169]:
x = np.arange(-100,100,.1)
plt.figure()
plt.plot(x, triangle(x))
plt.show()
In [170]:
type(np.array([1]))
Out[170]:
In [174]:
integrate.quad(triangle, -np.Inf, np.Inf)
Out[174]:
In [178]:
vec = metropolis_hastings(lambda f: 1/100*triangle(x), nsamples=10000)
In [181]:
x = np.arange(0,100,.1)
y = 1/100 * triangle(x)
plt.subplot(211)
plt.title('Metropolis-Hastings')
plt.plot(vec)
plt.subplot(212)
plt.hist(vec, bins=30,normed=1)
plt.plot(x,y,'ro')
plt.ylabel('Frequency')
plt.xlabel('x')
plt.xlim((-30,30))
plt.legend(('PDF','Samples'))
plt.show()
In [184]:
import pymc
In [185]:
@pymc.stochastic(dtype=int)
def switchpoint(value=1900, t_l=1851, t_h=1962):
"""The switchpoint for the rate of disaster occurrence."""
if value > t_h or value < t_l:
# Invalid values
return -np.inf
else:
# Uniform log-likelihood
return -np.log(t_h - t_l + 1)
In [189]:
from scipy.stats import rv_discrete
from scipy import stats
normdiscrete = stats.rv_discrete(values=(gridint, np.round(probs, decimals=7)), name='normdiscrete')
n_sample = 500
np.random.seed(87655678) # fix the seed for replicability
rvs = normdiscrete.rvs(size=n_sample)
rvsnd = rvs
f, l = np.histogram(rvs, bins=gridlimits)
sfreq = np.vstack([gridint, f, probs*n_sample]).T
print sfreq
In [219]:
from scipy import stats
class deterministic_gen(stats.rv_continuous):
def _pdf(self, x):
return 1/x
In [220]:
deterministic = deterministic_gen(name="deterministic")
In [222]:
deterministic.cdf(np.arange(1, 3, 0.5))
In [223]:
deterministic.pdf(np.arange(-3, 3, 0.5))
Out[223]:
In [224]:
deterministic.rvs(size=10)
In [195]:
deterministic?
In [239]:
def getDistribution(data):
kernel = stats.gaussian_kde(data)
class rv(stats.rv_continuous):
def _rvs(self, *x, **y):
# don't ask me why it's using self._size
# nor why I have to cast to int
return kernel.resample(int(self._size))
def _cdf(self, x):
return kernel.integrate_box_1d(-numpy.Inf, x)
def _pdf(self, x):
return kernel.evaluate(x)
return rv(name='kdedist')
In [240]:
f = getDistribution(np.random.randn(100))
In [243]:
f.rvs([1])
In [231]:
def getDistribution(data):
kernel = stats.gaussian_kde(data)
class rv(stats.rv_continuous):
def _cdf(self, x):
return kernel.integrate_box_1d(-numpy.Inf, x)
def _pdf(self, x):
return kernel.evaluate(x)
return rv(name='kdedist')
Out[231]:
In [244]:
kernel = stats.gaussian_kde(np.random.randn(100))
In [245]:
kernel
Out[245]:
In [246]:
type(kernel)
Out[246]:
In [247]:
kernel.resample(10)
Out[247]:
In [256]:
def getDistribution():
f = np.random.normal
class rv(stats.rv_continuous):
def _pdf(self, x):
return f(x)
return rv(name='kdedist')
In [257]:
a = getDistribution()
In [262]:
a.rvs(size=10)
In [264]:
import scipy.integrate as integrate
result = integrate.quad(lambda x: np.random.normal(x), 0, 4.5)